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Abstract 

This paper constitutes our initial effort in developing sparse grid discontinuous Galerkin 
(DG) methods for high-dimensional partial differential equations (PDEs). Over the past 
few decades, DG methods have gained popularity in many applications due to their dis¬ 
tinctive features. However, they are often deemed too costly because of the large number 
of degrees of freedom of the approximation space, which are the main bottleneck for simu¬ 
lations in high dimensions. In this paper, we develop sparse grid DG methods for elliptic 
equations with the aim of breaking the curse of dimensionality. Using a hierarchical 
basis representation, we construct a sparse finite element approximation space, reducing 
the degrees of freedom from the standard 0{h~^) to 0{h~^ \ log 2 for d-dimensional 

problems, where h is the uniform mesh size in each dimension. Our method, based on 
the interior penalty (IP) DG framework, can achieve accuracy of 0{h^\ log 2 h\^~^) in the 
energy norm, where k is the degree of polynomials used. Error estimates are provided 
and confirmed by numerical tests in multi-dimensions. 

Keywords: discontinuous Galerkin methods; interior penalty methods; sparse grid; high¬ 
dimensional partial differential equations. 

1 Introduction 

Elliptic equations have been used successfully in many areas of modeling, such as aeroacous- 
tics, electro-magnetism, oil recovery simulation, weather forecasting, etc. Their numerical 
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treatment includes the standard finite difference methods, collocation methods, Galerkin fi¬ 
nite element methods, least squares methods, among many others. Discontinuous Galerkin 
(DG) methods work well for purely hyperbolic problems due to the discontinuous nature 
of the solutions. Yet, they are also proven to be effective for elliptic equations. Since the 
1970s, the interior penalty (IP) DG methods have been proposed [HllMlEnilH] and generalized 
gaiiziES]. A review of DG methods for elliptic equations including discussions about other 
DG methods, such as the method of Bassi and Rebay |9] and the local DG (LDG) method 
[SQiiziiin], can be found in [0] . It is generally understood that the advantages of DG meth¬ 
ods lie in their flexibility in choosing the discretization meshes and the approximation spaces. 
However, the methods are often deemed too costly due to the large degrees of freedom of the 
approximation space. Such drawback is more prominent when dealing with high-dimensional 
equations arising from real-world applications, such as kinetic simulations, stochastic analysis, 
and mathematical modeling in finance or statistics. 

This paper concerns the computation of high-dimensional elliptic equations, for which the 
main challenge is commonly known as the curse of dimensionality P. This term refers to the 
fact that the computational cost and storage requirements scale as for a d-dimensional 

problem, where h denotes the mesh size in one coordinate direction. This challenge typically 
can not be resolved through barely increasing computational resources, and it requires the 
improvement of numerical techniques as well as efficient computational implementations. The 
sparse grid techniques [ISIEH], introduced by Zenger [60], have been developed as a major 
tool to break the curse of dimensionality of grid-based approaches. The idea relies on a 
tensor product hierarchical basis representation, which can reduce the degrees of freedom 
from 0{h~'^) to log 2 for d-dimensional problems without compromising much 

accuracy. The fundamentals of sparse grid techniques can be traced back to Smolyak [HT] for 
numerical integration, and they are closely related to hyperbolic cross [U |32] , boolean method 
[23] . discrete blending method [1^, and splitting extrapolation method m in the literature. 
The construction of the scheme is based on balancing the cost complexities and the accuracy 
of the scheme by seeking a proper truncation of the tensor product hierarchical bases, which 
can be formally derived by solving an optimization problem of cost/benefit ratios |3T]. Sparse 
grid techniques have been incorporated in collocation methods for high-dimensional stochastic 
differential equations [5H1 ETl HU 01], Galerkin finite element methods [60l Il5l 08], finite 
difference methods [301 [35] , finite volume methods 1371, and spectral methods [allEHl 09] EO] 
for high-dimensional PDEs. However, their potential has not yet been fully realized under the 
DG framework, which constitutes the main focus of this paper. 

Generally speaking, the design of DG methods has two essential ingredients: (1) a (weak) 
formulation of the methods based on the underlying PDE, (2) a suitable approximation space. 
The appropriate choice of the formulation of the scheme as well as approximation properties 
of the discrete space can guarantee the methods with desirable properties such as stability 
and convergence. In fact, DG formulations based on traditional piecewise polynomial space 
have been extensively studied for the last few decades, and the methods are rather mature 
for various applications. As for the approximation spaces, DG methods are more flexible 
compared with continuous finite element methods due to the lack of continuity requirement. 
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In recent years, the ideas of using nonstandard spaces such as non-polynomial spaces ISSllSl], 
or polynomial spaces with specihc properties such as locally divergence-free properties [20] 
have been explored. They are mainly driven by the needs to mimic particular properties of 
the exact solution. In this paper, we aim at incorporating the sparse grid approximation 
spaces in the DG framework to treat high-dimensional problems. As an initial effort in our 
investigation, we focus on high-dimensional elliptic problems and develop a sparse grid DG 
method for the following d-dimensional model problem 

- V • (KVn) = / in D = [0,1]", (1) 

with a Dirichlet boundary condition 

u = g on dfl, (2) 

where m(x) : G — )■ M is the unknown function and the matrix-valued coefficient function 
K = {ki,j)i<i,j<d is symmetric positive dehnite and bounded below and above uniformly, i.e. 
there exist two positive constant Kq,Ki, such that 

V X G G, A'o X ■ X < Kx ■ X < A'l X ■ X. 

Using a hierarchical basis representation, we construct a sparse hnite element approxima¬ 
tion space, reducing the degrees of freedom from the standard 0{h~'^) to 0{h~^ \ log 2 for 

d-dimensional problems, where h is the uniform mesh size in each dimension. Based on the 
IPDG framework and careful implementation, the new approximation space enables an effi¬ 
cient solver for the model equation. A priori error estimate shows that the proposed method 
is convergent in the order of 0{h^) up to the polylogarithmic term | log 2 in the energy 
norm when the solution is smooth enough, where k is the degree of the polynomial space. 
Although the method in this paper is restricted to problems on a box-shaped domain, it has 
the potential to be implemented on general domains by considering either a coordinate trans¬ 
formation or a more complex sparse hnite element space based on hierarchical decompositions 
on unstructured meshes. 

The rest of this paper is organized as follows: in Section we construct and analyze the 
DG approximation space on sparse grids. In Section we formulate the scheme and perform 
error analysis of the sparse grid IPDG method. Implementation details will be discussed to 
ensure computational efficiency. In Section numerical examples in multi-dimensions (up to 
d = 4) are provided to validate the accuracy and performance of the method. Gonclusions 
and future work are given in Section 

2 DG Finite Element Spaces on Sparse Grids 

In this section, we introduce the main ingredient of our algorithm: the DG hnite element 
space on sparse grids. We proceed in several steps. First, we review the standard piece- 
wise polynomial spaces in one dimension, and introduce a hierarchical decomposition with a 
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set of orthonormal bases. Second, we construct the sparse hnite element space from multi¬ 
dimensional multiwavelet bases as tensor products from the one-dimensional construction. 
Finally, we discuss some key features and approximation properties of the sparse hnite ele¬ 
ment space. Numerical tests are given to compare several dehnitions of sparse discontinuous 
hnite element spaces. 

2.1 Hierarchical Decomposition of Piecewise Polynomial Spaces in 
One Dimension 

In this subsection, we will introduce the hierarchical representation of piecewise polynomial 
spaces in one dimension. Many discussions here are closely related to the original Haar wavelet 
[36j, the multiwavelet bases constructed in laE], and multiresolution analysis (MRA) with 
DG methods for adaptivity [161IH EHl [381 [26] and constructions of trouble-cell indicator [53] 
for hyperbolic conservation laws. 

Without loss of generality, consider the interval fl = [0,1], we dehne the n-th level grid 
consisting of = (2“”j, 2“”(j -f 1)], j = 0,..., 2” — 1 with uniform cell size h = 2~'^. On this 
grid, we use || ■ to denote the broken Sobolev norm, i.e. -i = l^j=Q IrI ll-a/rjyj 

where is the standard Sobolev norm on /^. Similarly, | • denotes the broken 

Sobolev semi-norm, i.e. 

We can dehne 

to be the usual piecewise polynomials of degree at most k on this grid. Clearly, has degrees 
of freedom 2'^{k -f 1), and we notice that they have the nested structure for diherent values of 
n, 

< C C 14^ C 1 / 3 ^ C ■ ■ ■ 

Due to this nested structure, we can dehne the multiwavelet subspace n = 1, 2,... as the 
orthogonal complement of in iW with respect to the inner product on D, i.e. 

* n—1 ^ * n ^ n ^ * * n ^ n—1 

For notational convenience, we also denote the base space Wq ■= Vq, which consists of all 
polynomials of up to degree A; on [0,1]. The dimension of IF^ is 2”“^ (A; -|- 1) when n > 1, and 
A; -|- 1 when n = 0. Now, we have found a hierarchical representation of the standard piecewise 
polynomial space on grid {I^,j = 0,..., 2"' — 1} as = 0 o<n< 7 v ^n- We remark that the 
space Wll, for n > 1, represents the hner level details when the mesh becomes rehned and this 
is the key to the reduction in degrees of freedom in higher dimensions. 

For implementation purpose, we need to introduce basis functions for space W^. The case 
of n = 0 is trivial. It suffices to use a standard polynomial basis on D. For example, by 
using the scaled Legendre polynomials, we can easily obtain a set of orthonormal bases in Wq. 
When n> 1, we will use the orthonormal bases introduced in [2]. In particular, the bases for 
are dehned as 

hi{x) = 2^/‘^fi{2x-1), i = l,...,k + l 
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where {fi{x), i = 1,..., fc + 1} are functions supported on [—1,1] depending on k, with the 
following properties. 

1. The restriction of /* to the interval (0,1) is a polynomial of degree k. 

2. The function /j is extended to (—1,0) as an even or odd function according to the parity 
of i + k: 


f,(x) = (-l)-+*/,(-x). 
3. The bases have vanishing moments: 


fj{x)x^dx = 0 i = 0,1, ■ ■ ■ , j + k — 1. 


’-1 


4. The bases have the orthogonality and normality properties: 


fi{x)fj{x) dx =< fi, fj >= 6ij, i,j = l,...,k. 


'-1 


Functions {/j} can be computed by a repeated Gram-Schmidt algorithm. The particular 
form of {fi} up to fc = 4 are listed in Table 1 in [2]. Multiwavelet bases in [2] retain the 
orthonormal properties of wavelet bases for different hierarchical levels by introducing the 
basis for W^, n > 1 as 

= s'”-'’'" /i,(2”-U - j). » = 1,..., fc + 1. j = 0,.... 2”-‘ - 1. 

To make the notations consistent and compact, when n = 0, the bases v^q^x), i = 1,... k+1 
are dehned as the rescaled Legendre polynomials on [0,1] with orthonormal property. In 
summary, we have EE], 


where da/, 6nn', djj' are the Kronecker delta symbols. 

Next, we will discuss about the projection operators that are essential in finite element 
analysis. We dehne as the standard univariate projection operator from T^[0,1] to Vf. 
From there, we can introduce the increment projector, 

ifn>l 
Po^ ifu = 0, 




then we can see that 

W^ = QlL\0,l). 
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In summary, we arrive at the following identities for the hierarchical decomposition of the 
piecewise polynomial space and the projection operator: 

= 0 = E 

0<n<N 0<n<N 

The properties of naturally rely on P^. For the projection operators, we recall the 
following approximation results. 

Property 2.1 (Convergence Property of the Projection Operator [18j) For a func¬ 
tion V G 1), we have the convergence property of the projection as follows: for 

any integer t with 1 <t < min{p, k}, and s <t 1, 

\\PnV — ^ ^^||'y||ri'‘+i(o), (3) 

where Ck^s,t is a constant that depends on k,s,t, but not on n. 

From this property, using basic algebra, we can deduce that for n > 1, 

with 

= Cfc,,,t(l + 2-(*+^-'*)). (4) 

Finally, we want to remark on a subtle issue concerning the hierarchical decomposition and 
the multiwavelet space. If we use a projection other than P^, say is another projection 
from L^(f2) to Vfj. Then the increment projector ^ Q^. The space = (5^T^(0,1) is 
different from llP, and likewise for the multiwavelet bases. However, the hierarchical structure 
© IF ^ = Vjf still holds. In our discussion in the next subsection, we will highlight the 
implication of this statement and show that using different dehnitions of the projector and the 
increment space will not affect the dehnition of the sparse discontinuous hnite element space 
in multi-dimensions. 

2.2 Sparse Discontinuous Finite Element Spaces in Multi-dimensions 

In this subsection, we introduce the sparse hnite element space constructed from multi¬ 
dimensional multiwavelet bases obtained from tensor products of the one-dimensional bases 
in Section o 

For a d-dimensional problem, we consider the domain x = (xi,... ,Xd) G H = [0,1]“^. To 
facilitate the discussion, we hrst introduce some notations on the norms and operations of 
multi-indices in Ng, where No denotes the set of nonnegative integers. For a = (ai,..., ad) G 
Nq, we dehne the and norms 


d 



m=l 


ckIoo := max 

l<m<d 
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the component-wise arithmetic operations 

a ■/3 := (tti/Si,... c ■ a := (cai,..., ca^), 2“ := (2“i,..., 2"'*), 

and the relational operators 

cx < (3 am < /3m, 
cx < (3 OL < (3 and cx ^ f3. 

Now, for a multi-index \ = [li,. ■ ■ ,ld) ^ which indicates the level of the mesh in a multi¬ 
variate sense, we consider the standard rectangular grid hli with mesh size hi := (2“^h • • • > 2“^'*). 
We dehne the smallest size among all dimensions to be h = min{2“*b • • • > 2“^'*}, an elementary 
cell ^ = {x : Xm ^ (2"'™jm,2“'’"(jm -F 1)}, and 

Vf = {v : v(x) e 0<j<2'-l}, 

where P^(/j) denotes polynomials of degree up to k in each dimension on cell Jj. The space Vj' 
contains the traditional tensor-product piecewise polynomials used in the DG discretizations. 
Moreover, if we use equal rehnement of size 2~^ in each direction, we denote the space to be 
V^, and it consists of {2^{k -|- l))*^ degrees of freedom. 

The foundation of sparse grids is to use the tensor product of the one-dimensional hi¬ 
erarchical decomposition, and only chooses the most relevant bases to guarantee suitable 
approximation properties. To illustrate the main ideas, we introduce 


W? = W/" xWl" ■ ■ • X Wt 

*’^1 ll,Xi ^ ''12^X2 ^ ld,Xd'’ 


where Wk ^ corresponds to the space Wk dehned in the m-th dimension as dehned in the 
previous subsection. The number of degrees of freedom associated with Wj' is 

d 

dim{W\) = n dimiWt). 

m=l 

Based on the one-dimensional hierarchical decomposition, it is easy to see 

Vf = X ... X = 0 w!-, 


and 


Vt = VSx. X ■ ■ ■ X vs,., = 0 Wf, 


,<N 


The basis functions for Wj' can be dehned by a tensor product 


•= n jm = 0 ,..., max(0, 2 


Im 1 


l)i 1, . . . , /c ~h !• 


m=l 
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They form a set of orthonormal basis due to the property of the one-dimensional bases. 

Now, we are ready to introduce the sparse hnite element approximation space. In partic¬ 
ular, we dehne 


V jv • — 


0 wf. 


i<7V 


This dehnition is motivated by continuous hnite element space on sparse grid [U] . Instead of 
considering the standard piecewise polynomial space V^, which has exponential dependence 
on the dimension d, we shall use its subspace V^. This space has good approximation results 


(see Section 2.3 for details) with signihcantly reduced degrees of freedom. We remark that the 
key in the construction lies in the choice of condition for 1. Here we have taken it to be |l|i < N. 
As demonstrated in [T2] , other conditions are also possible to realize dimension reduction with 
aims of optimizing various types of approximation results. We defer the numerical study for 


the comparison of several different spaces to Section [2H 

Remark 2.2 We mentioned in the previous subsection that there are different ways to con¬ 
struct the increment space Wf. However, the choice will not affect the definition ofY^. This 


is because 


0 wf = 0 vf 


i<N 


1<N 


© wf, 

\Mi<N 


for any space Wf constructed from one dimensional increment space 
= V^. 


satisfying Vf_i © 


The next lemma will give a count of dimensions for the space V^. 
Lemma 2.3 The dimension of\^ is given by 

d—m—1 

dimfV’fi) = {k + (- 1 )'^-"* + 2^+™-'^+^ 



n=0 


N' 


n 


- 2 ) 


d—m—l—n 


+ 1 


Suppose there is an upper bound on the dimension d < do, then there exist constants Cdo,Cdo 
depending only on do, such that 


dr%N i\Td—l 


Cdfik + < dim{y%) < Cdfik + l)'^2^''iV' 

Proof: Due to the distinctiveness of the zero-th level in a mesh, we need to distinguish the 
zero-th level from other levels in the proof. Therefore, for each multi-index 1, we dehne the 
set of the dimensions with mesh level 0 as 1° = {z| /* = 0, z = 1,..., d} and the number of such 
dimensions as 

|l|o := # in 

Then 


dim{y%) = dim{ 0 Wf) 

|l|i<V 

d 

= E 

m=0 |l|o=m,|l|i<Ar 
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We will discuss two cases based on the value of |l|o. 

Case 1. 0 < |l|o = m < d — 1, i.e., there are m dimensions of level zero from the d 
dimensions of the multi-index 1. Clearly, 


dim{ 0 Wf) 

|l|o=m,|l|i<Ar 



dim{ 0 Wf). 

lt=0, l<t<m; /t>0, t>m 
|l|i<7V 


Since there are d — m dimensions that have meshes of level no less than one, we always have 
\\\\ > d — m. Dehne = (0, ■ • ■ , 0,1, • ■ ■ ,1) to be a vector in Nq whose hrst m dimensions 
are 0, and the rest d — m dimensions are 1. Then, by Lemma 3.6 in [T5] . 


dim{ 0 Wf) ={k + ly 


lt=0, /t>0, t>m 

Wi<N 


.(fc + l)" ^ 2 


lt=0, l<t<m; lt>0, t>m 
|l|i<Af 

N 

n—(d—m) 


E 


n=d—m 

N 


^t=0, /t>0, t>m 

|l|l=n 


.(fc + l)'' ^ 2 


n—(d—m) 


n=d—m 


n — 1 

{d — m) — 1 


= {k + 
= {k + 


N+m—d 


1)" E 


n=0 


l)d 



n + {d — m) — 1\ 
(d-m)-l ) 

d- 

I +m—d-^l 




Case 2. |l|o = d means 1 = 0. Therefore, 


dim{ 



Wf) = dim(W^) 


|l|o=rf, 1| 


i<Ar 


(fc + l)'". 


Finally, we combine both cases, and arrive at 


d 

dim{\%) = 0 Wf) 

m=0 |l|o=m,|l|i<V 


= (fc + l)'" 




•(- 2 ) 


d—m—l—n 
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We notice that 


d-l 

E 

m=0 
d-l 


m 


'\N-\-m—d-\-l 


d—m—1 

E 

n=0 


N' 


n 


- 2 ) 


d—m—l—n 




7n=0 






d 


m 


rn AJ—m 


2^N 


m=0 

= CdXN^~\^ + < CdXN‘^~\ 


where we use Cdg, Cd^ to denote generic constants that depend only on do, not on d, k, N. They 
may have different values in each occurrence throughout the proof. 

Similarly, 


d-l 


dim{\%) > dim{ 0 Wf) = {k + 1)'^ (-1)'^ + 2^-'^+^ ^ 

|l|o=0, l|i<V V 

>{k + > Cdoik + 


n=0 


(- 2 ) 


d—l—n 


In summary, there exist constants Cdfj,Cdo, such that 


Cd,{k + l)''2'^iV‘'-i < d*m(V^) < Cd,{k + iy2^N^-\ 


and we are done. ■ 

This lemma implies that, upon mesh refinement, the degrees of freedom for the sparse 
hnite element space will grow in the order of h~^ \ log 2 instead of h~‘^ for the traditional 
piecewise polynomial space. This translates into a significant reduction in computational cost 
when d is large. The next subsection will verify approximation property of the space V^. 


2.3 Approximation Results 

Given any one-dimensional projection operator ^ and its increment operator we can 
naturally define a projection from to as. 


Pn 


E e 

|l|i<Ar 


k 

—h,xi 


- ^di^d 


(5) 


where Qf denotes the operator Qf in the m-th dimension. In [18], continuous sparse 
hnite element space is used with streamline diffusion method to solve transport-dominated 
diffusion problem. The finite element space considered is := n for k > 1. 

The approximation properties for projector defined in (|^ in L‘^{kl), norms are 
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obtained when is a nnivariate projector onto nC°([0,1]) satisfying ([^. Because is 
a subset of V^, we can directly use the results in |18] and obtain an estimate of the projection 
error for P^. 

We introduce some notations about the seminorm of a function’s mixed derivatives. For 
any set I = {ii,... v} C {1, .. .d} with |/|=r, we dehne P to be the complement set of I in 
{1,... d}. For non-negative integers a, fd and set I, we dehne the seminorm 




E ■■■ E E E 

0<ai<Q 0<Or<a0</3i</3 0<Pd-k<l3 





Q!^d—r 


dx‘ 


Pd — r 
jd — r 


L2(n) 


and 


’\yt+i(n) ■= max max ( max Inlm+i.s.ifQ'i I . 
^ ^ s&{0,l} l<r<d \\I\=r ^ ’J 


We recall the following results in |1H] about the approximation properties in H^{Q) 

and H‘^{Q]s[) norm. 

Theorem 2.4 Let defined in (|^ be constructed from Pff which is a univariate projector 
onto Vfi nC°([0,1]) satisfying (|^, then for k > 1, any 1 < t < min{p, k}, there exist constant 
Cf, Ko(fc, t, N), Ki(/c, t, N) > 0, such that for any v G N > 1, d > 2, we have 


iP^v - v\ 




= \PmV - V 






for s = 0,1 norm and seminorm), and 


|PjvW - v\ 
where 


- 1 [d^/^ + d‘^Kfik,fiNY-^2-^) 


Cf.j. [d^^^Ki{k,t, NY + d‘^Ko{k,t, NY ^2 '^) 2 k>2 


k = 1, 


Kfik, t, iV) = I + Cyo,o, s - 0 

t 2Ckfl^t + Ckflfl, s — 1) 

and Cfc,o,i, Cfc,o,o o^re defined in ([^ and Q. Moreover, for s = 0, if 

~ 1) + Ckfifl < 1, 


there exists a positive constant Ct^k, such that K{k,t, s, N) < 1 for all N > l,d > 2 with 
N P 1 < Ct^kid — 1). 


This theorem implies a convergence rate of 0{h^^Y up to the polylogarithmic term | log 2 h\'^~^ 
in the norm, which is comparable with the traditional full grid approach when the function 
V retains enough smoothness (we actually do require higher regularity of the function v when 
compared with standard piecewise polynomial spaces). This estimate serves as the foundation 
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in the error estimates of IPDG method in the energy norm. We remark that when switching 
the projection operator (say projection or other projectors on to V^) the error bound may 


change. However, our estimates in Section 3.2 do not require any particular property from 
the projection. Therefore, it suffices for us to use this theorem directly from [IS]- In future 
work, we will establish convergence property of different projections onto V^, for /c > 0. This 
will be essential for error estimates for broader dehnitions of DG methods for other types of 
equations. 


2.4 Approximations from Various Sparse Discontinuous Finite El¬ 
ement Spaces: A Numerical Investigation 

In this subsection, we investigate the approximation properties of three sparse hnite element 
spaces from a numerical perspective. As one can imagine, there is an intricate relation between 
the number of bases used with the approximation properties of the spaces. Indeed, the choice 


adopted in Section [272] is not unique in sparse grids literature. Ideally, how to choose the “best” 
bases is an optimization problem of cost/beneht ratios [21], and the ultimate strategy may 
be an adaptive algorithm. Our preliminary study in this subsection compares three intuitive 
choices of spaces, including dehned in the previous subsections, and 


V'' — 

^ N ■ — 


0 wf, 


.|i<Ar+d-l 

|lloc<Ar 


which is a slightly larger space. contains the sparse continuous piecewise linear space i 


m 


when k > 1. Moreover, by a similar argument as in the proof of Lemma 2.3, its degrees 
of freedom, though slightly larger, are also of the order h~^ \ log 2 

The other space we consider retains not only sparsity with respect to h, but also with 
respect to the degree of polynomials (we call this p-sparsity). If we introduce 


X IT,^2 ... ^ ppfcd 

Xl l2,X2 '' 


with k = (fci. 


It is clear that 


wf = 0 

lkll<fc 

, fcrf), the new sparse space can be dehned as 


V jv • — 


0 wf. 


i<N 


ijk ^ ^ ^k 

^ N ^ ^ N ^ ^ N- 


A similar argument of Lemma 2.3 gives. 
Lemma 2.5 The dimension of\% is 

+ d\ f ^-4 / d 


dim{V%) = 


d 



d—m 


+ 2 


N+m—d-\-l 


d—m—1 

57 

n=0 


(- 2 ) 


d—m—l—n 


+ 1 
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Suppose there is an upper hound on the dimension d < do, then there exist constants Cdf,,Cdo 
depending only on do, such that 


Cdo 


k + d 




do 


k + d 
d 




Note that, when k and d are large, dim{\^) can be signihcantly smaller than dim{\^) and 
dim(V^). As for the implem enta tion of the space V^, we can no longer use the orthogonal 
hierarchical bases in Section 


2.1 


to represent functions in V^, since the vanishing-moment 
functions /j, which are used to construct the orthogonal hierarchical bases, have the same 
degrees. Instead, we can adopt another set of non-orthogonal hierarchical basis functions 
dehned as follows. We hrst construct the bases for W^. Again, the case of n = 0 is trivial. 
Similar to the orthogonal case, the bases for are dehned as 

hi{x) = 2^/‘^fi{2x-1), i = l,--- ,k + l, 

where fi, i = 1, ■ ■ ■ ,k + 1 are functions supported on [—1,1] and dehned by, 

fi{x)=x'-~^, xe(-l,0) 

and 

fi{x) = -x"~^, xe(o, 1). 

The bases for n>l are 

= 2^^-^^!^^ hi(T-^x - j), i = 1,..., fc + 1, j = 0,..., 2^-^ - 1. 


In summary, the basis functions in Wf (and V^) are 


: = 


n 

m=l 

|i|i<A: 


n-", (X. 

*771 ?*T71 


{xm), jm = 0,... ,max(0,2 


^771 1 


1 ). 


This set of basis is no longer orthonormal in the L?' sense, but a Gram-Schmidt algorithm can 
be adopted to achieve the purpose. 

Ideally, we should be able to use detailed analysis to compare the approximation results 
from the three spaces. However, for simplicity of implementation, here we only provide a 
numerical comparison for the approximation of a smooth function 


m(x) = exp ( Xi j , X G [0,1]'', 


771 = 1 


by the standard projections: P^, P^, and Pj( onto V^, V^, and V^, respectively. We 
measure various norms of the projection errors and the degrees of freedom for each space to 
hnd the best balance between accuracy and efficiency among the three choices. 
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In Tables 1 2 we report the projection errors ei = P^u — u, 62 = P^u — u and 63 = P^u — u 
in L“, and norms and the associated orders of accuracy for k = 2, d = 2,3. The 

detailed computational procedures as well as errors for other values of d and k are reported 
in [55]. For all three spaces, the degrees of freedom are all signihcantly less than {k + 
which is the degree of freedom for the full grid approximation. For space V^, [k + l)-th 
order of accuracy is clearly observed for the and errors, and k-th order of accuracy is 
observed for the error, while slight reduction of accuracy is observed for the L°° error. For 
space V^, k-th order of accuracy is observed for the error and slight reduction of accuracy 
is observed for and errors. However, we do observe about half-order to one order 
reduction of accuracy for the error. We remark that the magnitude of errors computed 


by space is larger than that by space for a fixed N. The apparent reason is that the 
degrees of freedom of are smaller than V^. However, when the degrees of freedom are 
comparable for the two spaces, the magnitude of errors computed by space is smaller than 

V^. Therefore, is preferred if the computational efficiency is concerned. For V^, which 
has the least degrees of freedom among the three sparse spaces, severe reduction of accuracy 
is observed for all error norms. Moreover, the magnitude of errors is significantly larger than 


that by space or with comparable degrees of freedom. 

In summary, based on the discussion, we will adopt space in the computation for 
elliptic equations in the next section. A further reduction in computational cost is possible by 
finding the optimal subset of by adaptive algorithms. We do not pursue this direction in 
this paper, and leave it to future study. 


3 IPDG Method on Sparse Grids for Elliptic Equations 

In this section, we solve the second order linear elliptic boundary value problem Q by IPDG 
methods with the sparse hnite element space V^. We will hrst formulate the scheme and 
discuss its implementation details, and then perform a priori error analysis. 

3.1 Formulation of the Scheme 

First, we introduce some basic notations about jumps and averages for piecewise functions 
defined on a grid VLn- Let F := UtsOjv union of the boundaries for all the elements 

in Hat and ^(F) ;= be the set of functions dehned on F. For any q G ^(F) 

and q G [5'(F)]'^, we define their averages {q'},{q} and jumps [g], [q] on the interior edges as 
follows. Suppose e is an interior edge shared by elements T+ and T_, we define the unit normal 
vectors n~^ and n~ on e pointing exterior to T+ and T_, and then 

[q] = g-n- +g+n+, {g} = i(g-+ g+), 

[q] = q-.n-+q+-n+, {q} = ^(q~ + q”^). 
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Table 1: Projection errors and orders of accuracy. SGDOF denotes the degrees of freedom 
of the sparse approximation space in use. FGDOF denotes the degrees of freedom of the full 
grid approximation space and is equal to 2^‘^(/c + 1)'^. d = 2. k = 2. 


N 

SGDOF 

FGDOF 

error 

order 

error 

order 

LOO 

order 

Lf i g];ror 

order 

^ N 

2 

108 

144 

2.82E-05 


4.36E-05 


4.05E-04 


2.26E-03 


3 

432 

576 

3.50E-06 

3.01 

5.47E-06 

2.99 

5.94E-05 

2.77 

5.66E-04 

2.00 

4 

720 

2304 

4.39E-07 

3.00 

6.86E-07 

2.99 

8.54E-06 

2.80 

1.42E-04 

2.00 

5 

1728 

9216 

5.51E-08 

3.00 

8.61E-08 

3.00 

1.21E-06 

2.81 

3.54E-05 

2.00 

6 

4032 

36864 

6.90E-09 

3.00 

1.08E-08 

3.00 

1.71E-07 

2.83 

8.84E-06 

2.00 

ijk 

''Af 

2 

72 

144 

3.51E-05 


5.23E-05 


6.48E-04 


2.41E-03 


3 

180 

576 

4.84E-06 

2.86 

7.26E-06 

2.85 

1.23E-04 

2.40 

6.08E-04 

1.99 

4 

432 

2304 

6.56E-07 

2.88 

9.96E-07 

2.87 

2.21E-05 

2.47 

1.53E-04 

1.99 

5 

1008 

9216 

8.81E-08 

2.90 

1.35E-07 

2.88 

3.77E-06 

2.55 

3.82E-05 

2.00 

6 

2304 

36864 

1.17E-08 

2.91 

1.81E-08 

2.90 

6.13E-07 

2.62 

9.55E-06 

2.00 

N 

2 

48 

144 

1.73E-03 


2.51E-03 


2.78E-02 


5.62E-02 


3 

120 

576 

4.86E-04 

1.83 

7.22E-04 

1.80 

1.09E-02 

1.35 

2.80E-02 

1.01 

4 

288 

2304 

1.35E-04 

1.85 

2.02E-04 

1.84 

3.99E-03 

1.46 

1.39E-02 

1.01 

5 

672 

9216 

3.68E-05 

1.88 

5.53E-05 

1.87 

1.37E-03 

1.55 

6.95E-03 

1.00 

6 

1536 

36864 

9.90E-06 

1.89 

1.49E-05 

1.89 

4.45E-04 

1.62 

3.47E-03 

1.00 
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Table 2: Projection errors and orders of accuracy. SGDOF denotes the degrees of freedom 
of the sparse approximation space in use. FGDOF denotes the degrees of freedom of the full 
grid approximation space and is equal to 2^‘^(/c + 1)'^. d = 3. k = 2. 


N 

SGDOF 

FGDOF 

error 

order 

error 

order 

LOO 

order 

Lfi grror 

order 

^ N 

2 

1188 

1728 

8.65E-06 


1.88E-05 


5.64E-04 


9.83E-04 


3 

4104 

13824 

1.08E-06 

3.00 

2.36E-06 

3.00 

8.11E-05 

2.80 

2.45E-04 

2.00 

4 

12096 

110592 

1.35E-07 

3.00 

2.95E-07 

3.00 

1.15E-05 

2.82 

6.12E-05 

2.00 

5 

32832 

884736 

1.69E-08 

3.00 

3.69E-08 

3.00 

1.65E-06 

2.80 

1.53E-05 

2.00 

6 

84672 

7077888 

2.12E-09 

3.00 

4.62E-09 

3.00 

2.40E-07 

2.78 

3.83E-06 

2.00 

ijk 

'' Af 

2 

351 

1728 

1.41E-05 


2.58E-05 


1.33E-03 


l.lOE-03 


3 

1026 

13824 

2.12E-06 

2.73 

3.86E-06 

2.74 

3.16E-04 

2.08 

2.80E-04 

1.98 

4 

2808 

110592 

3.15E-07 

2.75 

5.76E-07 

2.74 

7.07E-05 

2.16 

7.07E-05 

1.99 

5 

7344 

884736 

4.62E-08 

2.77 

8.56E-08 

2.75 

1.50E-05 

2.24 

1.77E-05 

1.99 

6 

18576 

7077888 

6.66E-09 

2.79 

1.26E-08 

2.76 

3.01E-06 

2.31 

4.44E-06 

2.00 

N 

2 

130 

1728 

4.50E-03 


6.59E-03 


1.13E-01 


9.32E-02 


3 

380 

13824 

1.82E-03 

1.30 

2.71E-03 

1.28 

7.36E-02 

0.62 

5.93E-02 

0.65 

4 

1040 

110592 

7.51E-04 

1.28 

l.llE-03 

1.29 

4.30E-02 

0.78 

3.91E-02 

0.60 

5 

2720 

884736 

3.10E-04 

1.28 

4.53E-04 

1.29 

2.32E-02 

0.89 

2.70E-02 

0.54 

6 

6880 

7077888 

1.30E-04 

1.25 

1.88E-04 

1.26 

1.18E-02 

0.97 

1.95E-02 

0.47 
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If e is a boundary edge, then we let 


[q] = qn, {q} = q, 


where n is the outward unit normal. 

Now we are ready to formulate the IPDG scheme for Q. We look for Uh G V^, such that 


B{u^,v) = L{v), VkeVJ, 


( 6 ) 


where 


and 


B{w, v) 



[u] ds 


^ [ {KVu} -Hds + J^T [N • H ds, 

eer de ggr ^ 


Hv) = /" /udx - / ( 
Jo, Jan ^ 


a 


KVu ■ n + —v) gds, 


(7) 

( 8 ) 


where a is a positive penalty parameter, h = 2 ^ is the uniform mesh size in each dimension. 
The IP methods developed in [5B1 [S] are popular methods for elliptic equations, but usually 
with standard piecewise polynomial space. Here, by using the sparse hnite element space, 
we can achieve a signihcant reduction in the size of the linear algebraic system (|^ especially 
when d is large. 

Several issues still have to be addressed to ensure real computational gains in implemen¬ 
tation. First of all, it is obvious that we need efficient algorithms to evaluate L{v) when v is 
taken to be the basis functions uli(x) in V^. A standard integration scheme will incur com¬ 
putational complexity with exponential dependence on d. To avoid this, in our algorithm, we 
take full advantage of numerical integration on sparse grids developed in the literature (SUET]. 
In particular, for each basis Ujj(x), its domain of dependence consists of 2^^ smooth patches. 
Therefore, when evaluating fvf j dx, we first divide this integral into 2'^ parts accordingly, 
and then for each patch, we implement a sparse grid integration with Gauss quadrature, and 
the number of levels employed in this calculation is taken as {Iquad — |l|i/2), where I quad 
is a fixed integer chosen large enough to guarantee sufficient accuracy. The evaluation of 

(KVu ■ n -h g ds can be performed in a similar fashion. 

Another important factor when we assemble the linear system is the so-called unidirec¬ 
tional principle [13]. To demonstrate the main ideas, we can see that evaluating f^</)(x)dx 
with 0(x) = (/>i(xi)... (/>d(xd) is equivalent to multiplication of one-dimensional integrals 
/[o 1 ] 01 dxi ■ • • /[Q (/>d dxd- Therefore, computing fvf j dx when f(x) is separable, i.e. f(x) = 
fi(xi)... fd(xd) or when f(x) is a sum of separable functions is straightforward, because we 
only need some small overhead to compute one-dimensional integrals and assemble them to 
obtain the multi-dimensional integrations. 
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The same discussion holds true for the computation of the bilinear term B{uh,v). For ex¬ 
ample, if we use a direct method, we need to evaluate B{w, v) for {w, v) being all possible basis 
functions in V^. From the dehnition of B{w,v), each matrix element will involve four multi¬ 
dimensional integrations. If K is separable (in particular, when K is a constant function), 
due to the unidirectional principle, the matrices can be assembled fast. When K is a general 
function, we need to compute true high dimensional integrals. This difficulty is identihed as 
one of the main challenging tasks for computing PDFs on sparse grids, see e.g. [HI [1]. In 
this case, we can either use the sparse grid integration procedure mentioned above or by a 
computational procedure outlined as follows. Assume K to be a smooth function, then we can 
hnd K/i = as the projection of K onto the sparse hnite element space and use 

K/i in place of K in the scheme. Notice that K.h is a sum of separable functions, therefore the 
computation of the bilinear term is accelerated as the unidirectional case. The reason we use a 
higher order sparse hnite element for projection of K.^ is to obtain exact evaluation of the vol¬ 
ume integral KVu ■ Vv dx. However, this process does change the values of two other terms 
up to approximation error of K^, — K. Another aspect we did not explore is the efficient solver 
for the linear algebraic system. In the literature, iterative methods have been proposed based 
on the semi-coarsening approach and its sparse grid extensions |l2l 031 |29l [Ml |33| 061112]. We 
leave those interesting implementation aspects to future study. 


3.2 Error Analysis 

This section contains error estimates of the IPDG method on sparse grids. Following [5|, we 
dehne the energy norm of a function v G by 


^ f iVnpdx + '^h 

tgUn eer 




We review some basic properties of the bilinear operator B{-,-). 

Lemma 3.1 (Orthogonality) Let u be the exact solution to ([^, and Uh he the numerical 
solution to (|^, then 

B{u-Uh,v)=0, Vn e V^. 

Proof. Using integration by parts, we can easily show B{u,v) = 0, Vn G V^. The Galerkin 
orthogonality immediately follows. ■ 


Lemma 3.2 (Boundedness [[51 IB]) There exists a positive constant Cf,, depending only on 
A'l, a, such that 

B{w,v) < G;,|||tc||| ■ |||n|||, G 

Lemma 3.3 (Stability [[51 [5]) When a is taken large enough, there exists a positive constant 
Cs depending only on Kq, such that 

B{v,v) >Cs\\\v\\\‘^, VugV^. 
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Theorem 3.4 Let u be the exact solution to Q, and Uh he the numerical solution to (|^. For 
k > 1, u e 1 < t < min{p, k}, N > 1, d > 2, we have 


u-Uh S 




where Cd is a constant that depends on d linearly. Cb^Cg are defined in Lemmas 3.2 and 3.3 


C* = max 


\ld^nl^-^ + + 2d^nl'^-^2-^^, \l + 2ci + 2d^4^-^2-^A , 


where Kg = Kg{k,t,N), s = 0,1 and are defined in Theorem 2.4 


Proof. Choose any function ui G V^, then we decompose the error into e = u — Uh = 


{u — uj) + {uj — Uh)- By Cea’s lemma, using Lemma 3.1, 3.2 and 3.3 


Cs|||m/ — w/iIlP < Bh{ui — Uh,ui — Uh) — Bh{ui — u,ui — Uh) < CfelHu — n/||| ■ HIm/ — Uh\ 


Therefore, 11 |m/ — m/i|| < ^||In — u/l11, and 


|e||| < |||n/ - M/i||| + |||m - M/Ill < ( 1 + ^ 1 lll"^ ~ 


We have 


Cb 


Cg 


e|||< l + TT inf |||m - M/Ill < 1 + — |||m - P^ulll, 


s J uj&V 


Cb 


Cg 


where the projection operator has been specihed in Theorem 


2.4 


Next, we need to bound the energy norm oi t] := u — P^u. Recall that 

■'cO.,ggp -J e K J .1 e 


Teftp 


The hrst term in the summation is 


eer 


[ \^v\^d^=\v\l^np 

TeOjv dT 


To bound the second and third terms, we use the trace inequalities |5]: 


II'^IIl2(9t) — Cd y ^ll0llL2('r) + j , V0 G H^{T) 


dfi 


dn 


L’^i&T) 
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where Cd is a generic constant that depends on d linearly. It may have different values in each 
occurrence in the proof. Sum over all the elements, we get 

eer ^ ^ 

f Ivf ds < Cd • 

eer ' “'e V / 

In summary, 


- {j^\\v\\L^{nN) + IvIh^Qn) + d 


By Theorem 2.4 


< Cdcl^ ^ + Sd^nf^ + ^2 2 fc > 2 

< Cdcl-i- + d^nj'^ + 2d + 2d‘^K,l'^~‘^2~‘^’^) 2“^^*|M|^t+i(j^), A; = 1 

where we have used the shorthand notation Ks = Ks{k, t,N), s = 0,1. Let’s define 

C* = max (\Jd‘^K^Q^~‘^ + Sd^nf^ + 2 #Kq‘’*“^ 2 “ 2 A^^ ^ d^Kg'^"^ + d'^nf^ + 2d + 2 #/?, 


2d-2r)-2N 
0 ^ 


Then, 

Therefore, we have proved the error estimate in the energy norm 


kill < '\/^Cfc,tC'*2 ^*|u|^t+i(n). 


This theorem implies a convergence rate of 0{h^) up to the polylogarithmic term | log 2 
in the energy norm when u is smooth enough. However, we do require more regularity of u 
compared with IPDG methods using traditional piecewise polynomial space. 

Remark 3.5 Proving convergence in norm with the standard duality argument will en¬ 
counter some difficulties in this framework. For example, let ip be the solution to the adjoint 
problem 

-V • (KVv^) = u — Uh on Q, p = 0 on dfl. 

Since G is convex, we have p G and || 95 ||r/ 2 (Q) < C\\u — Uh\\L^{n)- From the adjoint 

consistency of the IP method, we get 

lk-«/*llL2(o) = B{u-Uh,p). 

However, to proceed from here, we will need to define an interpolant of p: pi G and 

bound 11 |(p — (p/| 11. From our previous discussion, this will reguire a bound in | |(p| This 

is a stronger norm than the classical norm, and cannot be bounded by ||m — rLh\W 2 [Qy 
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4 Numerical Results 


In this section, we provide multi-dimensional numerical results to demonstrate the performance 
of our sparse grid IPDG scheme. 

4.1 Two-dimensional Results 

In this subsection, we gather computational results for two-dimensional case. The penalty 
constant in this subsection is chosen to be cr = 10 for k = 1, and cr = 20 for k = 2. 

Example 4.1 We solve the following two-dimensional problem with constant coefficient on 
G = [ 0 , 1 ] 2 . 

f-AM = 0, xgG, 

[u = sin(7rxi) X G dQ, 

for which the exact solution is u = sin( 7 ra:i) sePeme with k = 1 and k = 2 

on different levels of meshes. Numerical errors and orders of accuracy measured in L°° 

and norms are listed in Table We observe fc-th order of accuracy for norm, close to 
{k + l)-th order accuracy for norms and slightly less than {k + l)-th order accuracy for 

the L°° norm. The results agree with the error estimates performed in the previous section. 
The slight order reduction for and norms is similar to the one observed in Section 

1231 

In addition, we provide the sparsity patterns and the condition number of the stiffness 
matrices for k = 1 and fc = 2 in Figure [T] and Table From Table the number of nonzero 
elements scales as 0{SGD0F^'^), where SGDOF is the number of degrees of freedom of the 
space used. This is a denser matrix than the one generated from traditional piecewise polyno¬ 
mial space, for which one element can only interact with itself and its immediate neighbors. 
For sparse grids, the basis functions in are no longer locally defined due to the hierarchical 
structure. However, hxing any point x in the domain the number of bases with nonzero 
values at X is a constant for each level. In this sense, the interaction between the bases is still 
limited when compared with a global approximation scheme. 


Example 4.2 We solve the following two-dimensional problem with smooth variable coeffi¬ 
cient on H = [ 0 , 1 ]^. 

f-V ■ ((sin(a;ia; 2 ) + 1) Vn) =/, x G H, 

= 0, X G dQ. 

/ is a given function such that the exact solution is m = sin( 7 ra;i) sin( 7 ra; 2 ). As discussed in 
previous Section, when assembling the stiffness matrix, we hrst obtain the projection of 
coefficient K in space V^, then make use of the unidirectional principle to save computational 


21 





Table 3: Numerical errors and orders of accuracy for Example 
with fc = 1,2 and indicated N. 


4.1 


computed by the space 


N 

L} error 

order 

error 

order 

L°° error 

order 

error 

order 

k = 1 

3 

4.49E-03 


6.97E-03 


3.26E-02 


1.77E-01 


4 

1.18E-03 

1.93 

1.93E-03 

1.85 

9.71E-03 

1.75 

8.80E-02 

1.01 

5 

3.03E-04 

1.96 

5.09E-04 

1.92 

3.19E-03 

1.60 

4.36E-02 

1.01 

6 

7.68E-05 

1.98 

1.32E-04 

1.98 

9.68E-04 

1.72 

2.16E-02 

1.01 

k = 2 

3 

9.52E-05 


1.33E-04 


5.74E-04 


7.61E-03 


4 

1.42E-05 

2.75 

2.03E-05 

2.71 

9.65E-05 

2.57 

1.91E-03 

1.99 

5 

2.05E-06 

2.79 

3.02E-06 

2.75 

1.59E-05 

2.60 

4.78E-04 

2.00 

6 

2.89E-07 

2.83 

4.36E-07 

2.79 

2.66E-06 

2.58 

1.19E-04 

2.00 


Table 4: Sparsity and condition number of the stiffness matrix. Example 4T computed by 
the space with k = 1,2. SGDOF is the number of degrees of freedom used for the 
sparse grid DG scheme. NNZ is the number of nonzero elements in the stiffness matrix. 
Order=log(NNZ)/ log(SGDOF). 


N 

SGDOF 

NNZ 

Order 

Gondition Number 

k = 1 

3 

80 

992 

1.57 

3.58E+02 

4 

192 

3216 

1.54 

1.43E+03 

5 

448 

9168 

1.49 

5.68E+03 

6 

1024 

24144 

1.45 

2.26E+04 

k = 2 

3 

180 

3456 

1.57 

1.40E+03 

4 

432 

11124 

1.54 

5.49E+03 

5 

1008 

31596 

1.50 

2.16E+04 

6 

2304 

83028 

1.46 

8.58E+04 
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The sparsity patterns of the matrices computed by the space V| with 
fc = 1,2 and N = 6. Each dot represents a non-zero element in the stiffness matrix, (a): 
fc = 1 , (b): k = 2. 


Figure 1: Example 


4.1 


cost. The numerical results are provided in Table The conclusion is very similar to the 
constant coefficient case Example 4T We observe fc-th order of accuracy for norm, close 
to {k + l)-th order of accuracy for L^,L^ norms and slightly less than {k + l)-th order of 
accuracy for the norm. 


Example 4.3 We solve the following two-dimensional problem with discontinuous coefficient 
on = [ 0 , 1 ]^. 


I-V-((sign((xi- 0 . 5 )(a; 2 - 0 . 5 )) + 2)VM) = /, x G fl, 

= 0, X e dfi. 

/ is a given function such that the exact solution is u = sin( 7 ra:i) sin( 7 rx 2 ). In this example, 
the coefficient has four constant panels on but the solution itself is smooth. We want to use 
this example to demonstrate the performance of the proposed scheme when solving problems 
with discontinuous coefficient. We provide the numerical results with k = 1 and k = 2 in Table 
1^ Notice that in this example, the discontinuity of the coefficient is along the cell boundaries 
for all our computational grids, and the convergence rates are similar to the previous two 
examples. 


4.2 Three-dimensional Results 


In this subsection, we gather computational results for three-dimensional elliptic equations. 
The penalty constant in this subsection is chosen to be cr = 15 for k = 1, and cr = 30 for 


k = 2. 


23 

































Table 5: Numerical errors and orders of accuracy for Example 
with fc = 1,2 and indicated N. 


4.2 


computed by the space 


N 

error 

order 

L?' error 

order 

L°° error 

order 

error 

order 

k = 1 

3 

1.30E-02 


1.65E-03 


4.50E-02 


3.37E-01 


4 

3.18E-03 

2.03 

4.08E-03 

2.01 

1.34E-02 

1.74 

1.66E-01 

1.02 

5 

7.81E-04 

2.02 

l.OlE-03 

2.01 

4.43E-03 

1.60 

8.26E-02 

1.01 

6 

1.94E-04 

2.01 

2.55E-04 

1.99 

1.48E-03 

1.58 

4.11E-02 

1.00 

k = 2 

3 

1.77E-04 


2.17E-04 


5.74E-04 


1.35E-02 


4 

2.71E-05 

2.70 

3.37E-05 

2.69 

l.OlE-04 

2.51 

3.37E-03 

2.00 

5 

3.99E-06 

2.76 

5.08E-06 

2.73 

1.91E-05 

2.40 

8.41E-04 

2.00 

6 

5.67E-07 

2.82 

7.37E-07 

2.78 

2.99E-06 

2.67 

2.10E-04 

2.00 


Table 6: 
with k = 


Numerical errors and orders of accuracy for Example 
1,2 and indicated N. 


4.3 


computed by the space 


N 

L} error 

order 

error 

order 

L°° error 

order 

error 

order 

k = 1 

3 

1.24E-02 


1.57E-02 


4.55E-02 


3.33E-01 


4 

3.07E-03 

2.02 

3.94E-03 

2.00 

1.36E-02 

1.75 

1.66E-01 

1.00 

5 

7.58E-04 

2.02 

9.78E-04 

2.01 

4.50E-03 

1.59 

8.32E-02 

1.00 

6 

1.89E-04 

2.00 

2.46E-04 

1.99 

1.50E-03 

1.58 

4.16E-02 

1.00 

k = 2 

3 

1.96E-04 


2.59E-04 


1.21E-03 


1.56E-02 


4 

2.72E-05 

2.85 

3.50E-05 

2.89 

1.55E-04 

2.96 

3.70E-03 

2.08 

5 

3.85E-06 

2.82 

4.94E-06 

2.82 

2.05E-05 

2.92 

8.93E-04 

2.05 

6 

5.36E-07 

2.84 

7.02E-07 

2.82 

3.01E-06 

2.82 

2.19E-04 

2.03 
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Table 7: Numerical errors and orders of accuracy for Example 
with k = 1,2 and indicated N. 


4.4 


computed by the space 


N 

error 

order 

error 

order 

L°° error 

order 

error 

order 

k = 1 

3 

1.29E-02 


2.19E-02 


1.09E-01 


2.85E-01 


4 

4.05E-03 

1.67 

6.98E-03 

1.65 

4.75E-02 

1.20 

1.44E-01 

0.98 

5 

1.07E-03 

1.92 

1.94E-03 

1.85 

2.34E-02 

1.02 

7.02E-02 

1.04 

6 

2.76E-04 

1.96 

5.22E-04 

1.89 

8.44E-03 

1.47 

3.39E-02 

1.05 

k = 2 

3 

1.41E-04 


2.06E-04 


1.26E-03 


1.05E-02 


4 

2.51E-05 

2.49 

3.80E-05 

2.44 

3.35E-04 

1.91 

2.72E-03 

1.95 

5 

4.18E-06 

2.59 

6.49E-06 

2.55 

6.51E-05 

2.36 

6.87E-04 

1.98 

6 

6.69E-07 

2.64 

1.06E-06 

2.62 

1.09E-05 

2.58 

1.72E-04 

2.00 


Example 4.4 We solve the following three-dimensional problem with constant coefficient on 

n = [0,1]T 


—An = 0, X G 

ti = sm(7rx,) x e dSl, 


( 12 ) 


where the exact solution is n = sin(7ra;i) sin(7ra:2) We provide the numerical results 

with k = 1 and A; = 2 in Table For this three-dimensional example, we obtain fc-th order 
for the norm, and close to {k + l)-th order for the and norm, but L°° error seems to 
degrade to {k + |)-th order. However, upon closer examination, similar behaviors have been 

for the projection error of a smooth function onto V^. The sparsity 


2.4 


observed in Section 

patterns of the stiffness matrices for k = 1 and k = 2 are reported in Figure and Table 
From Table we can see the stiffness matrices scale less than the two-dimensional examples, 
i.e., the number of nonzero elements scales as 0{SGD0F^ '^), where SGDOF is the number 
of degrees of freedom of the space used. 


Example 4.5 We solve the following three-dimensional problem with smooth variable coeffi¬ 
cient on H = [0,1]^. 


-V ■ ((sin(a;ia;2X3) + 1) Vu) = /, 

u = 0, 


X G H, 

X G d^l, 


(13) 


/ is a given function such that the exact solution is u = sin(7rxi) sin(7ra;2) sin(7rx3). Again, 
for such a non-separable coefficient K, the projection onto space is obtained and the 
unidirectional principle is applied when we assemble the stiffness matrix. We provide the 
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Table 8: Sparsity and condition nnmber of the stiffness matrix. Example 4A compnted by 
the space with k = 1,2. SGDOF is the nnmber of degrees of freedom nsed for the 
sparse grid DG scheme. NNZ is the nnmber of nonzero elements in the stiffness matrix. 
Order=log(NNZ)/ log(SGDOF). 


N 

SGDOF 

NNZ 

Order 

Gondition Nnmber 

k = 1 

3 

304 

3760 

1.43 

3.73E+02 

4 

832 

14080 

1.42 

1.51E+03 

5 

2176 

45760 

1.39 

5.97E+03 

6 

5504 

135872 

1.37 

2.36E+04 

k = 2 

3 

1026 

20250 

1.43 

1.58E+03 

4 

2808 

74628 

1.41 

5.98E+03 

5 

7344 

240516 

1.39 

2.32E+04 

6 

18576 

710532 

1.37 

9.15E+04 


(a) 




0 6000 12000 18000 


(b) 


Fignre 2: 


Example 


4.4 


(a) /c = 1 and (b) k = 


, The sparsity patterns of the matrices compnted by the space Vg with 
2. Each dot represents a non-zero element in the stiffness matrix. 
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Table 9: Numerical errors and orders of accuracy for Example 
with fc = 1,2 and indicated N. 


4.5 


computed by the space 


N 

error 

order 

error 

order 

error 

order 

error 

order 

k = 1 

3 

2.64E-02 


3.40E-02 


1.55E-01 


4.32E-01 


4 

6.23E-03 

2.08 

8.58E-03 

1.98 

3.54E-02 

2.13 

2.04E-01 

1.08 

5 

1.49E-03 

2.06 

2.10E-03 

2.03 

2.07E-02 

0.77 

9.82E-02 

1.06 

6 

3.68E-04 

2.02 

5.32E-04 

1.98 

7.58E-03 

1.45 

4.80E-02 

1.03 

k = 2 

3 

1.63E-04 


2.05E-04 


8.24E-04 


1.19E-02 


4 

2.88E-05 

2.50 

3.66E-05 

2.48 

1.63E-04 

2.34 

3.00E-03 

1.98 

5 

4.72E-06 

2.61 

6.06E-06 

2.60 

2.73E-05 

2.58 

7.54E-04 

2.00 

6 

7.42E-07 

2.67 

9.58E-07 

2.66 

5.80E-06 

2.23 

1.88E-04 

2.00 


numerical results for k = 1 and A; = 2 in Table The 
including the comparable convergence rate and similar 
errors. 


conclusion of Example AA holds here, 
order reduction of the and 
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4.3 Four-dimensional Results 

In this subsection, we gather computational results for four-dimensional case. The penalty 
constant in this subsection is chosen to be cr = 30 for k = 1, and cr = 60 for k = 2. 

Example 4.6 We solve the following four-dimensional problem with constant variable coef¬ 
ficient on = [0,1]^. 


-Am = 0 


u = sin(7rxi) sin(7rx2) sin(7rx3) 


X G ff, 
X G dkl 


(14) 


where the exact solution is m = sin(7rxi) sin(7rx2) sin(7rx3) 

We list numerical errors and orders in Table For this four-dimensional example, we 
obtain fc-th order for the norm, and close to {k -|- |)-th order for the and norms. 
The L°° error seems to fluctuate. However, we expect the order of accuracy in L°° norm to 
grow upon refinement. The numerical results when fixing two of the coordinates are plotted 
in Figure]^ for k = 1 and k = 2 when N = 4. The left and right figures are obtained with 
degrees of freedom equal to 3072 and 15552, respectively. We can see the higher order scheme 
offers a much more accurate solution than the lower order scheme, and this is also verified 
The sparsity patterns and condition number of the stiffness matrices for k = 1 
are reported in Figure and Table IT 


by Table [10 
and k = 2 


From Table 11 


we can see the stiffness 
matrices scale less than the two-dimensional and three-dimensional examples, i.e., the number 
of nonzero elements scales like 0{SGDOF^'^^), where SGDOF is the degree of freedom of 
the space used. 



Figure 3: Example 4^ (a) fc = 1 ; (b) A; = 2. iV = 4. Plotted along xi = 0.4930, X 2 = 0 
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Table 10: 
with k 


Numerical errors and orders of accuracy for Example 4.6 computed by the space 
= 1,2 and indicated N. 


N 

error 

order 

error 

order 

LOO gJ-J-QJ. 

order 

error 

order 

k = 1 

3 

2.44E-02 


4.22E-02 


3.31E-01 


3.91E-01 


4 

1.08E-02 

1.18 

2.08E-02 

1.02 

1.16E-01 

1.51 

2.37E-01 

0.73 

5 

3.68E-03 

1.54 

7.15E-03 

1.54 

9.33E-02 

0.31 

1.22E-01 

0.96 

k = 2 

2 

8.21E-04 


1.34E-03 


l.llE-02 


4.20E-02 


3 

1.76E-04 

2.22 

2.79E-04 

2.27 

2.76E-03 

2.00 

1.20E-02 

1.81 

4 

3.32E-05 

2.40 

5.39E-05 

2.37 

8.76E-04 

1.66 

3.18E-03 

1.91 


Table 11: Sparsity and condition number of the stiffness matrix. Example 4.6 computed 
by the space with k = 1,2. SGDOF is the number of degrees of freedom used for the 
sparse grid DG scheme. NNZ is the number of nonzero elements in the stiffness matrix. 
Order=log(NNZ)/ log(SGDOF). 


N 

SGDOF 

NNZ 

Order 

Gondition Number 

k = 1 

3 

1008 

12272 

1.36 

4.27E+02 

4 

3072 

51712 

1.35 

2.26E+03 

5 

8832 

187008 

1.34 

9.27E+03 




k = 2 


2 

1539 

19683 

1.35 

7.40E+02 

3 

5103 

102303 

1.35 

2.62E+03 

4 

15552 

420336 

1.34 

9.72E+03 




(a) 


(b) 


Figure 4: Example 4.6 The sparsity patterns of the matrices computed by the space V 4 with 


{a) k = 1 and (b) k = 2. Each dot represents a non-zero element in the stiffness matrix. 
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Example 4.7 We solve the following four-dimensional problem with smooth variable coeffi¬ 
cient on = [0,1]^. 


-V ■ ((sin(a;ia;2a;3a;4) -1 - 1) Vm) = /, 

M = 0, 


X E f2, 
X E dQ, 


(15) 


/ is a given function such that the exact solution is m = sin(7ra;i) sin(7rx2) sin(7ra;3) sin(7ra;4). 
We provide the numerical results with k = 1 and k = 2 in Table The conclusion is similar 


to Example 4.6 


Table 12: 
with k 


Numerical errors and orders of accuracy for Example 4T computed by the space 
= 1,2 and indicated N. 


N 

error 

order 

error 

order 

LOO 

order 

Ri gj-j-or 

order 

k = 1 

3 

6.15E-02 


8.97E-02 


2.94E-01 


6.67E-01 


4 

1.89E-02 

1.70 

2.63E-02 

1.77 

2.54E-01 

0.21 

3.20E-01 

1.06 

5 

4.51E-03 

2.07 

6.80E-03 

1.95 

7.15E-02 

1.83 

1.45E-01 

1.14 

k = 2 

2 

8.38E-04 


1.09E-03 


3.49E-03 


3.74E-02 


3 

1.62E-04 

2.37 

2.13E-04 

2.36 

1.34E-03 

1.38 

l.OlE-02 

1.90 

4 

2.97E-05 

2.44 

3.91E-05 

2.45 

3.80E-04 

1.82 

2.57E-03 

1.97 
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5 Conclusion 


In this paper, we develop a sparse grid IPDG method for elliptic equations on a box shaped 
domain. The method features a sparse hnite element space which scales as 0{h~^ \ log2 h\^~^) 
for d-dimensional problems, translating into a signihcant cost reduction when d is large. The 
traditional IPDG formulation can be readily used, and when combined with approximation 
results for the new space, equip the scheme with error estimate that is only slightly deteriorated 
in the energy norm. Numerical results in up to four dimensions are shown to validate the 
performance of the method. Future work includes the study of adaptive algorithm for less 
regular solutions, extension of the method to more general domains, and developing sparse 
grid DG methods for other types of high-dimensional equations. 
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